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ABSTRACT 


A two-dimensional finite difference model was developed to cal- 
culate the temperature distribution during the fusion welding process. 
The model is based upon the enthalpy formulation for phase change 
processes. Changes in thermal properties during the process have 
been accounted for. Temperature variations with respect to time have 
been calculated for specific locations in the cross-section of the work 
piece away from the edges for specified boundary conditions. These 
computations allow the estimation of cooling rates and the extent of 
the fusion zone. The model could be used in real-time monitoring of 


the welding process. 
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I. INTRODUCTION 


A. BACKGROUND 

The welding process in most applications is largely dependent 
upon the skill of the welder with regard to the quality and accuracy of 
the weld. Control of the voltage and current settings and the speed at 
which the torch advances can significantly affect the welding process. 
An experienced welder can adjust these parameters based upon the 
conditions he observes during the process. As the experience level 
varies from welder to welder, the quality of the weld can vary also. An 
unsatisfactory weld will result in expensive and time-consuming 
rework, assuming that the flaw is detected and identified. If the flaw 
remains undetected, catastrophic failure may occur at some future 
time. Neither of these situations is desirable. If a reliable automated 
system can be developed that produces uniform and repeatable results 
for a specified set of conditions, the reliability of the welding process 
will increase and the long-term costs will decline. 

To support automating the process, a simple means of monitoring 
the weld is desired. One way of sensing what is occurring in the weld 
process is to monitor the heat flow patterns on the material being 
welded. If the heat flow patterns for a high-quality weld are estab- 
lished, the amount of deviation of the actual weld heat flow patterns 
from the established values can be used as a measure of the quality of 


the weld. When the quality level reaches a level considered 


unsatisfactory, a control system can be used to provide real-time cor- 


rective action. This area has been the subject of several studies. 


B. PREVIOUS STUDIES 

The first significant studies of heat flow in welding were con- 
ducted by Rosenthal starting in the late 1930s. Rosenthal investigated 
both point and line heat sources as idealizations of the welding pro- 
cess. Results of such analyses are only applicable far away from the 
heat sources. Between 1940 and 1980, there were several studies of 
heat transfer during the welding process, as described by Koichi 
Masubuchi [Ref. 1]. 

W. H. Giedt explored the problem of modeling a control system 
for the Gas-Tungsten-Arc welding process [Ref. 2]. A two-dimensional 
model originally proposed in Reference 3 is assumed to be a reason- 
able approximation of the heat flow in a work piece being welded. The 
two-dimensional model was chosen over the more-exact three-dimen- 
sional model, primarily due to the reduced computational time that is 
required. For such a system to be practical, the welding system should 
be portable. A portable computer implies reduced computational 
speed from a large mainframe computer, thus the emphasis on 
reduced computational time. 

Goldak, Chakravarti, and Bibby have taken the welding model a 
step further by modifying the grid pattern used [Ref. 4]. In order to 
improve the ability of the model to accommodate both shallow and 
deep penetration welds, the authors developed a non-axisymmetric, 


three-dimensional, heat-source model. To accommodate the new heat 


source, the grid pattern changed from rectangular to ellipsoid. This 
approach resulted in values that were closer to experimental values for 
both deep and shallow penetration welds than previous models. How- 
ever, the model became significantly more complicated to achieve the 
increased accuracy. 

Kou and Le present a three-dimensional heat flow model for a 
pipe undergoing a seam weld [Ref. 5]. This model addressed the prob- 
lem of weld pool convection by using an “effective” thermal conduc- 
tivity to account for the effect of convection within the weld pool on 
the heat flow. 

Oreper and Szekely delve more fully into the circulation patterns 
that occur within a weld pool [Ref. 6]. Their analysis centers around 
three key components. The heat flow must be modeled in both solid 
and molten regions. In the molten regions, fluid flow must take into 
account buoyancy, surface tension, and electromagnetic forces. 

The predictions from a computational model must be experimen- 
tally verified. Masubuchi [Ref. 1] describes several studies on tempera- 
ture measurement using thermocouples. Khan, et al., discuss the use 
of an infrared sensor to monitor weld surface temperatures [Ref. 7]. 
Infrared imaging of the surface area around the point of application of 
the weld torch provided a means of monitoring temperatures 
throughout this small area. As pointed out by them, the resolution of 
such a system must be sufficient to discriminate between disturbances 
that would affect weld quality and system noise. Additionally, the time 


constant of the system must be sufficiently small to be able to react to 


a detected problem before the torch has passed the region of the flaw. 
Unless the problem is corrected while the area is still molten, it will 
not be correctable. The need for a realistic model of the welding pro- 
cess for real-time control is clear from the study by Khan, et al. 


[Ref. 7]. 


C. MODELING THE HEAT FLOW 

In order for a real-time control system to be effective, the model 
must provide accurate data about the heat flow in the work piece. The 
model developed in the present study calculates enthalpy and 
temperature values in a vertical cross-sectional plane of the work 
piece. The data is delivered as temperature at each nodal location. 
Surface temperature data are obtained from a time history of such 
cross-sections. 

The flow of heat within the work piece includes several complex 
interactions. In the solid phase, conduction, surface convection, and 
radiation constitute the heat transfer nodes. In the molten region, liq- 
uid phase material is subject to added effects caused by weld pool cir- 
culation [Ref. 6]. The need for simplifying assumptions becomes 
evident when computational speed is necessary, as in real-time con- 


trol applications. These are discussed in the following chapter. 


II. TWO-DIMENSIONAL FINITE DIFFERENCE MODEL 


A. OBJECTIVES 

The present study examines a two-dimensional model which 
simulates the heat flow in a cross-section of a work piece during the 
autogenous (no filler metal) welding of a butt joint with no contact 
resistance (see Figure 1). The model starts with all temperatures in 
the cross-section at ambient conditions. As the leading edge of the 
welding torch reaches the plane of the cross-section, energy is added 
to the cross-section from the top surface, initially at an increasing 
rate. As the center of the torch passes the cross-section, the energy 
added decreases until none is being added when the trailing edge of 
the torch has passed the cross-section. The entire surface area of the 
slab is considered to have free convection with the ambient environ- 
ment. Additionally, the entire surface area is able to radiate to the 
environment. 

The model takes into consideration the effects of phase changes 
within the work piece. This is accomplished by using an enthalpy for- 
mulation. The model also accommodates various geometric dimen- 
sions, welding parameters (current, voltage, efficiency, torch speed of 
advance, etc.), and material properties. The accuracy to which the 
location of the reported temperature is determined could be varied by 


changing the number of nodes at which temperature is calculated. 
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Figure 1. Work Piece Schematic 


A two-dimensional finite difference model was chosen in this 


study due to its simplicity and faster computational speed. 


B. ENTHALPY FORMULATION 
Welding, by its very nature, involves a phase change from solid to 
liquid and eventually back to solid. Models considering uniform spe- 
cific heat have no way of accounting for the latent heat of fusion 
absorbed by the material as it melts once it has reached melting tem- 
perature. The temperature of the molten material will not rise above 
the melting temperature until its enthalpy has reached the sum of the 
melting enthalpy and the latent heat of fusion. By using the enthalpy at 
each node instead of the temperature in the storage term of the 
energy equation, phase changes can be easily accommodated and the 
temperature can be determined at each time step. This technique has 
been employed in the current computation. It is described in detail in 
Reierence 9), 
1. Derivation of Nodal Finite Difference Equations 
There are several categories of nodal finite difference equa- 
tions for the model of a rectangular slab. These include interior nodes, 
surface nodes not subject to a heat flux, surface nodes subjected to a 
heat flux, and surface corner nodes. The equations for the interior 
nodes, a surface node with heat flux, and a surface corner node are 
developed in detail next. The equations for all the categories can be 


found in Appendix A. 


All of the categories start with an appropriate energy balance 
over the control volume being considered. For an interior control vol- 


ume, this leads to: 

Of ot ee CR yale! 

ax(k Gx) * aur aioe Oy 
Since the thermal conductivity is a function of temperature, it will vary 
with location as the temperature at that location changes. This change 
in thermal conductivity with location is handled by calculating an 
interface conductivity for a position halfway between the node being 
investigated and each of its adjacent neighbors. This value is referred 
to as the harmonic mean value [Ref. 9]. For example, the conductivity 


between node (i,j) and node ({i-1,j) is found as follows: 
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}(AXL+AX) #(k,,AXL+k,_, AX) (3) 
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kk, (AXL + AX) 
i-1.j 


Sey led oe, ik (4) 


where AX is distance in X direction across node (i,j), and 


AXy is the distance in X direction from center of node (i,j) to 
center of node (i-1,j). 


A similar procedure is followed for each value of conductivity appro- 
priate to each of the four boundaries of the nodal square (see Figure 2). 
The equations for each conductivity can be found in Appendix A. The 
energy balance equation for an interior node becomes: 


reeGh = 1 J)- 10d) *AY  K2 « (Tle = 1) -T(LJ)) * AX 
AXL it AYT 


K3*(T(1+1.J) -TOJ)) AY | K4 *(TUd + 1)- TOW) #AX 
33 AXR AYB (5) 
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H™(1J) = At *(K1 *T(I- 1) * AY/AXL + KZ *T(LJ - 1) * AX/AYT 
+K3*T(I+1,.J) *AY/AXR +K4 *T(Ld + 1) * AX/AYB 
(6) 
~ (Kl * AY/AXL + K2 * AX/AYT +K3 + AY/AXR 


+K4 * AX/AYB) * T(I,J)) + H (J) 


The energy balance for a top wall node exposed to flux Q becomes (see 


Figure 3): 
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Figure 2. Control Volume Diagram for Interior Points 
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H**'(1,J) = At *(Q * AX +K1 * T(I- 1,J) * AY A2 * AXL) 
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The energy balance for a surface corner node becomes (see Figure 4): 
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Figure 4. Control Volume Diagram for NW Corner Point 
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C. SPECIAL CONSIDERATIONS 
1. Non-Uniform Grid Spacing 

The welding operation presents a situation where a highly 
concentrated heat source is applied to a small region of a surface at a 
much lower temperature. This situation results in large temperature 
gradients in the vicinity of the torch. In order to increase the accuracy 
of the temperature values at any given location, the concentration of 
nodes in the affected area has been increased. In the horizontal direc- 
tion, one-half of the nodes were located in the length directly under 
the heat source. The remaining half of the nodes were equally spaced 
over the distance from each side edge to the edge of the heat source. 
Similarly, one-half of the nodes were vertically positioned in the top 
20 percent of the thickness. 

The result is a fine grid mesh in the region directly under the 
heat source where the temperature changes from node to node are 
the greatest. As the distance from the heat source increases, the 
spacing between adjacent nodes increases, as seen in Figure 5. 

2. Modeling the Welding Torch Arc 

As seen in Figure 6, the heat source has been modeled as a 
Gaussian distribution. The rate at which energy is imparted from the 
torch to the surface is assumed to equal the voltage multiplied by the 


current multiplied by the welding efficiency. 
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Figure 5. Non-Uniform Grid Illustration 
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Figure 6. Welding Torch Arc 
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where n, the welding efficiency, is an empirical factor to account for 
the energy losses from the arc and the weld pool surface. 

Equation 11 also gives rise to a Gaussian distribution of the 
heat flux q”. The maximum value of q” occurs at the center line of the 
distribution. This peak flux value q% was established such that the 
integral of the flux values over the area the torch covers (a circle of 


diameter FW, the flame width) is equal to 7VI. 





CnVvl 
q., = 4 where C = mnt 
(=) (12) 
2 


The constant C was chosen such that the flux was reduced to a value of 
1 percent of the center line value by the time it had reached the 
maximum distance from the center line (FW/2). 

The model of the arc is designed to simulate the passage of 
the torch over a given cross-sectional plane. Equation 13 is the final 
equation taking both horizontal position and the passage of time into 


account. 


2 2 
” DQ a= - -( ) 
q’=q’,e C(L /2-XTOT ) e ClFWw/2 TIMTOT )(v) (13) 


The first exponential term takes into account the horizontal 
location. The upper left-hand corner is always the origin. XTOT is the 
total distance traveled in the x direction. L/2-XTOT represents the 


distance away from the center line. C is the same value as in equation 
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12. Far away from the center line, the argument of the exponential is a 
very large negative number, hence q” approaches zero. As XTOT reaches 
the edge of the flame, q” equals one percent of qg. As XTOT reaches 
L/2, q” = qo. The squaring of the distance allows the equation to decay 
as XTOT increases above L/2 in a mirror image of how it was built up. 

The second exponential term is very similar to the first. 
However, this term accounts for the approach and passage of the torch 
over the cross-section of interest. FW/2 is the distance from the cen- 
ter line of the torch to the edge of the flame. TIMTOT * V accounts for 
the distance that the torch has traveled from its initial starting point. 
The computations were started when the leading edge of the torch 
was 1/4 FW away from the cross-section of interest. Choosing this 
location allowed for observing the effects of the approach of the weld- 
ing torch. As was the case in the horizontal direction, the flux value 
builds up to a maximum when the torch is directly over the cross-sec- 
tion of interest and decays off in a mirror image of the way it built up. 

By multiplying the two exponential terms together, the true 
bell shape of the flux distribution is achieved. 


D. SOLUTION PROCEDURE 
This model of heat flow has been developed using the explicit 
technique. This implies that the enthalpy values for the new time level 
are calculated from known enthalpy values at the previous time level. 
The enthalpy at each node was initialized based upon its control 
volume and the ambient temperature (293K). The enthalpy values for 


each node at time level n + 1 are determined using the equations in 
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Appendix A. The temperature at each node at time level n + 1 can 
then be determined through equation 14. 


oe Hod) 


me) a Graxae a4) 


where _ pp = density of material 
AXAY = control volume term for 2D system 
Cp = specific heat 

The values for specific heat (Cp) are also temperature dependent. 
In order to use an explicit marching scheme, Cp was evaluated at tem- 
peratures calculated at time level n. Due to the small time steps used, 
the error due to this approximation is likely to be very small. 

Once the temperature at a node has reached the melting temper- 
ature, equation 6 must be modified. The latent heat of fusion of the 
material must be accounted for. In order to get a correct temperature 
value once the material is molten, the latent heat of fusion must be 
subtracted from the enthalpy of the control volume, as seen in the fol- 


lowing equations. 


ENFUS = (LATENT HEAT OF FUSION) pAXAY (15) 
ENMELT = (TMELT) CppAXAY (16) 


For H+1(1,J) > ENMELT + ENFUS 


n+l 


H (I,J) -ENFUS 


n+] 
T”* (I,J) = p C, AKAY (17) 


Wes: 


For ENMELT < H"+1(1.J) < ENMELT + ENFUS 
[n+l = TMELT (18) 


ro 


II. RESULTS 


A. PARAMETERS 

The model was run in tests of various combinations of the variable 
parameters. Test runs were made for two materials, steel and alu- 
minum. The specific properties of the materials can be found in 
Appendix B [Ref. 11]. 

Because the model used the explicit method for solving the finite 
difference equations, computational stability was a matter of concern. 
For a two-dimensional problem, the Fourier number for an interior 
node must always be less than or equal to 0.25 to maintain stability 
[Ref. 12] (see equation 19). Expressions similar to equation 19 can be 


derived from the boundary nodes. 


(AX) pC, (AX) (19) 


The result of equation 19 is to tie the precision of location of any given 
node, the AX value, to the time step required to maintain stability. As 
the precision increases (AX decreases), the time step decreases by the 
square of the change in AX. This situation puts a premium on the 
minimum required accuracy in order to maintain reasonable time 
steps. 

The parameters for the standard work pieces and torch settings 


were established as follows: 
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Width 0.35m Voltage 28 Volts 

Thickness 0.15m Efficiency 0.80 

Current 170 Amperes Speed of Advance 2.54 x 10°3 m/sec 
These dimensions and parameters were chosen to allow for compari- 
son with results obtained from an instrumented steel slab with similar 
parameters performed by Kihara, et al., and described in Reference 1. 

The first four computational runs for the two materials were made 

with varying spacing. The fifth run was made with an increased torch 
speed of advance. The sixth run was made with a reduced current. The 


seventh run was made with a smaller work piece. 


B. DISCUSSION 

The data from all runs can be found in Appendix D. 

The first four sets of results show the effect of grid spacing upon 
the computations. In all cases, as the grid spacing became smaller, the 
temperature values are likely to become more accurate. The 40 x 40 
grid spacing was taken as the optimum value since the highest tem- 
perature enter line value at time 8.8 seconds) was less than one 
percent different from the value obtained with a 50 x 50 grid spacing. 
The increase in accuracy obtained with the 50 x 50 grid spacing, 
though minimal, had a high cost in computational speed. The 50 x 50 
grid required approximately twice the amount of time to run as did 
the 40 x 40 grid. 

The data followed expected trends in that the steel heated up toa 
much higher temperature but the heat-affected zone was much 


smaller than for aluminum. This is appropriate because the thermal 
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diffusivity of aluminum (aa) =~ 97 x 106 m2/sec) is approximately five 
times greater than the thermal diffusivity of steel (are =~ 18 x 106 
m2/sec). 

With current reduced by 10 percent from 170 Amps to 153 Amps 
and using a 41 x 41 grid, the maximum temperature rise achieved on 
the center line of the cross-section was reduced by approximately 10 
percent in both the aluminum (11 percent) and steel (9.5 percent) 
samples. The depth of penetration of the molten region was also 
reduced as expected. Less energy being put into the samples resulted 
in less penetration and a lower peak temperature. 

The increase in speed of advance of the torch by 10 percent from 
2.54 x 10°3 m/sec to 2.794 x 10°3 m/sec resulted in a reduction of the 
increase in maximum center line temperature of approximately five 
percent. Again, the depth of penetration was reduced. In this case, the 
energy level was not reduced but it was not allowed to act on the 
cross-section as long. The net result was still less energy being input 
to the cross-section. 

As a check of the model against experimental data, the 41 x 41 
baseline grid was allowed to continue to run until 45 seconds had 
passed since the torch passed the cross-section (total time: 60.8 sec- 
onds). The center-line surface temperatures after the torch had 
passed were plotted against the data available from the instrumented 
slab experiments [Ref. 1]. The correlation between the model and the 


experimental data is good, especially since the material properties of 
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the experimental steel had to be approximated with data from the 


most similar alloy for which data were available (see Figure 7). 
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Figure 7. Comparison of Computational Results and Experimental 
Data From Kihara for Center-Line Surface Temperatures After the 
Heat Source Has Completely Passed the Center Line Surface Node 
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IV. CONCLUSIONS 


The model developed here is an adequate simulation of the heat 
flow during fusion welding at locations far away from the ends of the 
work piece. Improvement can be made by addressing the convective 
heat transfer within the molten weld pool. 

The model is based upon an explicit technique and hence is 
dependent upon the availability of a fast computer in order to handle 
small dimensions with good precision. To be practical for use on a 
mobile system, a fast, unconditionally stable model must be developed. 

The model should be expanded to three dimensions in order to 
be reliable in all locations of the work piece, not just in the interior 
regions. 

Actual experiments with instrumented work pieces should be 


performed to verify the data obtained from the model. 
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APPENDIX A 
MODEL EQUATIONS 


CONDUCTIVITY EQUATIONS 


k k,, *(AXR + AX) 


* 
i-1,j 


Kl = k CoC)... kk\ _oa_aoa—Xa—>;- 

1. iy ky, #AXR+k,,, *AX 
25a k, 4. *K,,* (AYT + AY) 

9. Beedle AV Tks co AY 
ee k,4,* Ky * (AXR + AX) 

5 ty ky * AXR +k, ,* AX 
" K, a * Ky, * (AYB + AY) 
K4= = 

an eos Kh SAY 


NODAL EQUATIONS 


EN(I,J): enthalpy at node (i,j) from previous iteration [J/m] 


EN1(I,J): enthalpy at node (i,j) for current iteration [J/m] 


TIME: time step [sec] 

H: convection (film) coefficient [W/m2 K] 
ae]: ambient temperature [K] 

Gg. flux [W/m?] 
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1. NW Corner 
ENI1(I,J) =TIME #((AX + AY) *H *TI/2+K3 *T(I+1,J) 
* AY (2 * AXR) + K4 * T(IJ + 1) * AXA2 * AYB) 
— ((AX + AY) *H/2+ K3 * AYA(2 * AXR) + K4 
+ AX (2 * AYB)) * T(1,J)) + EN(I,J) 
2. SW Corner 
ENI(1,J) = TIME * ((AX + AY) *H + TI/2 +K3 *T(1 + 1,J) 
* AY (2 * AXR) + K2 * T(1J - 1) * AXA2 * AYT) 
- ((AX + AY) *H/2+ K2 * AYA(2 * AXR) + K2 
* AX (2 * AYT)) * T(I,J)) + EN(I,J) 
3. Left Wall (Not Corners) 
ENI(1Jd) =TIME*(AY *H * T1+K3 *T(1+ 1,J) * AY/AXR 
+K2 *T(IJ—- 1) *AXA2 * AYT) + K4 *T(1,J +1) 
* AX (2 * AYB) -(H * AY + K3 * AY/AXR+ K2 


* AX [(2 * AYT) + K4 * AX /(2 * AYB)) * T(1J)) + EN(I,J) 


IAS 


4. NE Corner 
EN1(I,J) = TIME * ((AX + AY) *H * TI/2 + Kl * T(I - 1,J) 
* AY A2 * AXL) + K4*T(I,J +1) * AXA2 *AYB) 
— ((AX + AY) * H/2+ Kl * AYA2 * AXL) +K4 


* AX (2 * AYB)) * T(1,J)) + EN(I,J) 


5. Top Wall—Flux (Not Corners) 
ENI(I,J) = TIME *(Q * AX + Kl * T(I - 1,J) * AYA2 * AXL) 
+K3*T(0+1,J) *AYA2 *AXR) +K4 *T(LJ+)) 
* AX /AYB — (K1* AYA2 * AXL) + K3 * AYA2 * AXR) 


+K4 * AX/AYB) * T(I,J)) + EN(IJ) 


6. Top Wall—Flux/No Flux Interface 
ENL(1J) = TIME *(Q * AX/2 + AX *H *TI/2 + K1*T(I- 1,J) 
* AY (2 * AXL) + K3 * T(I+1,J) * AYA2 * AXR) 
+K4*T(I,J + 1) *AX/AYB -(H * AX/2+K1 
* AY (2 * AXL) + K3 * AYA2 * AXR) + K4 


+* AX /AYB) * T(1,J)) + EN(I,J) 
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7. Top Wall—No Flux 
ENI(1LJ) =TIME *(AX *H *TI+K1 *T(I- 1,J) 
+ AY (2 * AXL) + K3 * T(1+1,J) * AYA2 * AXR) 
+K4+*T(I,J +1) *AX/AYB -(AX *H+Kl 
# AY (2 * AXL) + K3 * AYA2 * AXR) + K4 


* AX /AYB) * T(1,J)) + EN(I,J) 


8. SE Corner 
ENI(1J) = TIME *((AX + AY) *H * TI/2+ Kl * T(1- 1,J) 
+ AY (2 * AXL) + K2 * T(J — 1) * AXA2 * AYT) 
— ((AX + AY) *H/2+ Kl * AYA2 * AXL) +K2 


* AX (2 * AYT)) * T(I,J)) + EN(I,J) 


9. Right Wall (Not Corners) 

EN1(1J) =TIME*(AY *H * TI+K1 * T(1—- 1,J) * AY /AXL 
+K2 *T(I,J -— 1) *AXA2 * AYT) + K4 *T(I,J + 1) 
* AX (2 * AYB) - (AY *H+K1* AY/AXL + K2 


* AX (2 * AYT) + K4 * AXf(2 * AYB)) * T(I,J)) + EN(I,J) 
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10. Bottom Wall (Not Corners) 

ENI(1I,J) = TIME * (AX *H *TI+K1 *T(I—- 1,J) * AYA2 * AXL) 
+K3*T(I+1,J) *AY/AXR +K2 *T(I,J - 1) * AX/AYT 
~ (AX *H +K1* AYA2 * AXL) + K3 * AY/(2 * AXR) 
+K2 * AX/AYT)*T(1J)) + EN(I,J) 

11. Interior Points 

ENI(1J) = TIME *(K1 *T(I—- 1,J) * AY/AXL + K3 * T(1+ 1,J) 

* AY /AXR + K2 * T(I,J — 1) * AX/AYT +K4 
*T(IJ +1) * AX/AYB — (K1 * AY /AXL +K3 * AY /AXR 


+K2 *AX/AYT +K4 * AX/AYB) * T(1,J)) + EN(I,J) 


ox) 


APPENDIX B 
MATERIAL DATA 


1. Aluminum 
99.996% pure 
Melting Temperature = 933.5K 
LHF = 3.98 x 105 J/kg 


a. Conductivity 
Temperature [K] Conductivity [W/m K] 
27 3:2 239 
373-2 234 
473.2 228 
Da.3 221 
673-2 215 
Ele, 208 
873.2 201 
973.2 103 
1073.2 12) 


Approximated by: 
K = 248.918398 - .02075127 T - 6.324 x 10°5 T2 + 2.78 x 10°8 TS 
[W/m K] 


b. Specific Heat 


Temperature [K] Cp [cal/g K] 
373.0 2.2832 x 10-1 
473.0 2.3795 x 10°! 
573.0 2.4759 x 10°! 
673.0 2.5760 x 10-1 
773.0 2.6723 x 1071 
873.0 2.7687 x 1071 
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Approximated by: 
Cp = 803.1932 + .4071 T [J/kg K] 


2. Steel 


Alloy: .23C 1.50 Mn .12 Si 
Melting Temperature = 1810 K 
LHF = 2.473 x 105 J/kg 


a. Conductivity 
Temperature [°K] Conductivity [W/m K] 
273.2 46.0 
323.2 46.0 
S732 46.0 
423.2 45.6 


Approximated by: 
K = 46.7358 - .0024 T [W/m K] 


b. Specific Heat 
Alloy: .13C .61 Mn .2 Si 


Temperature [K] Cp [cal/g K] 
273 .1039 
293 .1061 
323 .1109 
373 .1180 
423 1221 
473 .1260 
5Z3 1291 
avo else 
623 13S) 
6% .1429 
723 .1501 
CIES: -1582 
823 .1671 
873 .1802 
923 .1953 
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Temperature [K] Cp [cal/g K] 


973 Palen 
1013 2414 
1036 .2916 
1053 .2130 
1073 lshen 
1113 ebSi7 1 
1141 pez. 1 
aS .1570 
1223 .1561 
1273 .1570 


Approximated by: 
Cp = 418.458 + .2115 T [J/kg K] 


The numerical data was obtained from Reference 11. 
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APPENDIX C 
COMPUTER PROGRAM OF MODEL 


Pe 8s ote xf et <td fe os 3 ss oe ote ate 
ais alk iz ale he eke ois she dhe chs oh ois oe a oie ke Ss eh he oh ak ook oe 


Sedioh OGRA CALCLLATES THE TEMPERATURE 
* DISTRIBUTION IN A FLAT PLATE WHICH IS 

[we EOseED TO CONVECTIVE AND RADIATIVE 

* LOSSES ON ALL SIDES AND HAS A GAUSSIAN 

els tRIbL TION HEAT SOLRCE PASS DOWN THE 

eer TEREINE OF THE PLATE: 


IE Me StON 1(051.051).EN1(931,051),OU TPT(6,031), 

& TES 1(6),EN(050,051)79(051,U51).S(051,031), 

& KMELT(051.051) 

Ree ON DYER .LWANTVL.QO.AIBI.CI].D1,TIME.RHO.C, 
eee ee NOT TO ENIEN OUP, TESTS. 11K1,K2.KS, 
See LT TNF LT. PW.O0.V C2 KNIELT.EFF.VOLT.CURR. 
Peele ENT US. TINTOT,FACTOR DXNL-DNR.DY T.DYB, 
me NTOT.DALDS: 

Pave oie FINIShLCNIT.BAR.N.Z,PLEXEE LEX WFE.U, 

ae J1FEAG 

COnP Om il) )A1,BICiDI,T,.K1,K2.K3,K4,K MELT,DXL, 

& DXR.DYT,DYB 

OPEN(UNIT=1,FILE= OUTPY’) 

Al = 46.7358 

B1 = -.002400148 

C1l=9 

D0 


whe ah shen oh the aha ah th me th 2h ok OE BAR PARA METE RS 279 22 2% 2 2 He 2 ah ee Be oe ae ae oe oe Be ae He ae ae 
te, 


L=0.26250D0 
W =0.11250D0 
Err —( spe 
VOLT =028.D0 
CURR= 170.D0 
RHOQ= 73sb2. 03 
He oe 
FACTOR — 0 
Ne UO 
TIMORG = .02000D0 
TIMES TIXIORG 
PW. = 2.0008 -02 
"= 2.540L-03 
TWEE — Isto 


te gts Se ate he ate o's wo she oMe ate a ote ate ate gh gM a's a's ako gs p%s gto ake glo sto te ate ats sto te ss 5% gfe 5% 3! to sPs s% ale ole ate gle ote 3 sig ge s! yo gz ate x’ she oo o% oe o' 
ot she ee Se ole te ait gg cis ch she ae ote he ole he he ais ate ahs ate ahs ale he als oh ote he ote ate Se ate oh ahs ae oe ake oe oe oh oh Se 0 oe 0 Ooh ie eo 9 oe oe oo oe 


FLUXE= NINT(.5¢N + 1)-ENT(25®) 
FLUN2=NIN Bee te es 
QTOT=0 
WRITE(*,*) FLUXI,FLUX2 
C2=-LOG(.01) (FW 2)**2 
QO = EFF*VOLT*CURR®*C2 3.141593 
TI=293.D0 
DO 100 J=1,N+1 
DO 105 1=1,N+1 
KMELT(I,J)= 1.D0 
Saba XW CORNER 

IF ((J.EQ.1).AND.(1.EQ.1)) THEN 

S(I.J) = .25D0 

sees ener SW CORNER 

ELSE IF ((J.EQ.N+1).AND.(I.EQ.1)) THEN 
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Seelj= 25D0 


sseaeens LEFT WALL (NOT CORNERS) 
ELSE IF (1.EQ.1) THEN 
SUJ)=.5D0 
seseeees NE CORNER 


ELSE IF ((J.EQ.1).AND.(1LEQ.N+1)) THEN 
Sf) ==5D0 
seaecese TOP WALL (NOT CORNERS) 
ELSE IF ((J.EQ.1).AND.(ILGT.FLUX1).AND.(I.LT.FLUX2)) 
& THEN 
SUJ)=.SD0 
seseesin SE CORNER 
ELSE IF (J.EQ.N+ 1).AND(LEQ.N+1)) THEN 
SJ) =.28D0 
sasneaus RIGHT WALL (NOT CORNERS) 
eget O.\ +1) THEN 


Spl) = Ss D0 
aon BOTTOM WALL (NOT CORNERS) 
ese ir (J.EOQ.N+ 1) THEN 
sa) = -.SD0 
ia INTERIOR POINTS 
ELSE 
oy) — 120 
ESD'LF 


W1J)= 295.0 
IF (.GT.FLUKI1).AND(LLT.FLUX2)) THEN 
DXN = FW (FLUX2-FLUN1) 
DE1=DX 
Eee EO LUN!).OR(LEQIPLUX2)) THEN 
Dy —(@ (FE LUN2-FEUX1) + (L-FW) (N+ FLUXI-FLUX2)):2 
ELSE 
x =(L-PW) (\+ FLUX1I-FLUX2) 
DX2=DN 
END I 


105 
100 
109 


110 


[FLT NIN eaten 
DY =0.2*W NINT(.5*(N +1) 
ELSE IF (J.EQ.NINT(.5*(X+1))) THEN 
DY =(.2*W/NINT(.5°(N 4 1) +.8°W,(N-NINT(.5*(N 4 1))))'2 
ELSE 
DY =0.8*W(N-NINT(.5*(N 41) 
END IF 
EN(I,J)= 7.85E + 03*480.428*293*DX*DY*S(I.J) 
T9(IJ) = 293 
CONTINUE 
CONMINGE 
FORMAT(1LX,21F5.2) 
TIMTOT=0 
CNT=0 
CNT=CNT#1 
TIMTOT=TIMTOT+ TIME 
DO 200 J=1,N+1 
XTOT=0 
DO 205 [=1,N4+1 
ca — ie 
iF (GT FLU) AND (Eh Lie) aie 
DX =FW(FLUX2-FLUN1) 
IF (I-1.EQ.FLUX1} THEN 
DXL = ((L-FW)\(N + FLUX1-FLUX2)+ FW.(FLUX2-FLUX1))'2 
DXR=DX 
ELSE IF (1+ 1.EQ.FLUX2) THEN 
DXL=DX 
DXR=((L-FW) (N+ FLUXI-FLUX2)+ FW,(FLUX2-FLUX]))'2 
ELSE 
DXL=DX 
DXR=DX 
END IF 
ELSE IF (I.EQ.FLUX!) THEN 
DX =(FW(FLUX2-FLUX1)+(L-FW)'CX + FLUX1-FLUX2))2 
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Deer = (Loewen + FLUX }-PDUX?2) 
DNR =F WFLUN2-FLUNI) 
ELSEuieme@. FLU2)anHEX 
Di (F WALL M2-FLUX1)+(L-FW)'(N + FLUXI-FLUX2)) 2 
DAL = POAELUS2-FLUX1) 
eek =( 1  (N + FLUX1-FLUX2) 
ELSE 
x —(LoPWICNS FLUNFEUNQ) 
ibd Mee: FLUNI@HEN 
DXL=DX 
DNR=(FW (FLUX2-FLUNI)+(L-FW) (N+ FLUX1-FLUX2))'2 
ELSE IF (I-1.EQ.FLUX2) THEN 
rer (FWCFLUN2-FLUNI)+(L-FW)QN + FLUXE-FLUN?)) 2 
DXR=DX 
ELSE 
DXL=DX 
DNR=DX 
END IF 
EXDTE 
IF J.LT.NINT(.3*(N = 1))) THEN 
DY =0.2*W NINT(.3°(N +1) 
IF (J¥1.EQ.NINT(.5#(N + 1))) THEN 
DYT=DY 
DYB=(.S*W (N-NINT(.S4(N + 1) +.22W 
&  ININTCS(N S19) 2 
ELSE 
DYT=DY 
DYB=DY 
END IF 
ELSE IF (J.EQ.NINT(.5*(N + 1))) THEN 
DY =(.2*W NINT(.54(N+ 1))+.8?W (X-NINT(.5%(N + 1))) 2 
DYT=.2*W NINT(.52(N +1) 
DYB=.S*W (N-NINT(.52(N+1))) 
ELSE 


Bi 


DY =0.8*W.(N-NINT(.54(N + 1))) 
DYT=DY 
DYB=DY 
END IF 
I=] 
Jl=J 
K1=0 
K2=0 
K3=0 
K4=0 
H = HC + 4.0*5.67E-08*0.5*T(I,J)**3 
eae fete NW CORNER 
IF ((1.EQ.1)AND.I.EQ.1)) THEN 
NL COND: 
CAI GOn. D4 
ENKIJ)= TIME*((DX + DY)*H*TI 24+ K3*T(14 1)" 
& DY(2*DXR) +K4*T(LJ+1)*DX (2*DYB)-((DX + DY) 
& *H2+K3*DY'(2*DXR)+ K4*DX (2*DYB)) 
&) #4(1dj) + ENa oy 
sesh a SW CORNER 
ELSE IF ((J.EQ.N+1).AND.(I-EQ.1)) THEN 
CALECON DE 
CALE CONDS 
EN1IJ)= TIME*(DX + DY)*H*TI 2+ K3*T(1 + 1,J)* 
& DY(2*DXR)+ K2*T(1J-1)*DX (2*DYT)-((DX + DY) 
& *H 2+K3*DY:(2*DNR)+ K2*DX (2*DYT)) 
& y+ ENC LID 
fetes te LEFT WALL (NOT CORNERS) 
ELSE IF (1.EQ.1) THEN 
CAI COnD?2 
CALL GON. 
CALL COND4 
ENI(1,J)= TIME*(DY*H*TI + K3*T(1 + 1,J)*DY/DXR 
& +K2*T(LJ-1)*DX Q2*DYT)+K4*T(LJ+1) 
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& *DX (2*DYB)-(H*DY+K3*DY DXR+K2*DX (2*DYT) 
+ K4*DX (2*DYB))*T(1.))+ EN(IJ) 
aneenes NE CORNER 
ELSE IF ((J.EQ.1). AND.(ILEQ.N + 1)) THEN 
CALL COND 
CALE COND® 
ENI(1.J)= TIME*((DX + DY)*H*TI 2+ K1*T(I-1,J) 
& *DY (2*DNL)+K4*T(IJ+ 1)*DX (2*DYB)-((DX + DY) 
& *H2+KI1*DY \2*DNL)+K4*DX (2*DYB)) 
& *T(IJ))+EN(1J) 
anbaanee TOW r Lex (NOT CORNERS) 
Piette EQ AND (GT. FUUNI)AND(ULT.FLUX2)) 
& THEN 
IF (EXP(-C2*(L 2-XTOT)**2).LT.0.01) GO TO 999 
IF (EXP(-C2*(2.00*FW 2-TIMTOT#V)**2).LT.0.01) GO TO 
& 999 
CALL CONDI 
CALL COND3 
CANN CONDS 
Q= Q0*EXP(-C2°(L 2-XTOT)**2)*EXP(-C2*(2.00% FW 2 
& -TIMTOT#V)**2) 
ENI(13)= TIME*(Q*DX + K1®T(E-1,J)*DY (2*DXL) 
& +K3*T(1+1J)*DY (2°DXR)+ K4#T(1J+ 1)*DX DYB-(K1 
& *DY (2°DNL)+ K3*DY (2*?DXR)+K4*DX DYB) 
& *TUJ)y+ ENS) 
seeeeans TOP WALL-FLUX NO FLUX INTFC 
ELSE IF ((J.EQ.1). AND.((ILEQ.FLUX1).OR. 
& (LLEQ.FLUX2)),AND.(FLUXI.NE.FLUX2)) THEN 
IF (EXP(-C2*(L 2-XTOT)**2).LT.0.01) GO TO 999 
IF (EXP(-C2*(2.00*FW 2-TIMTOT#V)**2).LT.0.01) GO TO 
& 999 
@NUL, COND! 
CALL COND3 
CRUL COND4 
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Q= QO*EXP(-C2#(L/2-XTOT)**2)*EXP(-C2*(2.00* FW 2 
& -TIMTOT*V)*#2) 
ENI(1.J)= TIME*(Q#DX '2+ DX*H*TI 2+ KI1*T(I-1,J) 
&  *DY (2*DNL)+ K3*T(1+1J)*DY (2*DNR)+ K4*T(1,J +1) 
& DX 'DYB-(H*DX 2+ K1*DY,(2*DXL)+ K3*DY/(2*DXR) 
&  +K4*DX DYB)*T(I,J))+ ENGI) 
ELSE IF (J.£EQ.1),AND.((LEQ.FLUXI).OR. 
& (LEQ.FLUX2)).AND.(FLUX1.EQ.FLUX2)) THEN 
IF (EXP(-C2*(L!2-XTOT)*#*2).LT.0.01) GO TO 999 
IF (EXP(-C2*(2.00*F W 2-TIMTOT#*V)**2).LT.0.01) GO TO 
& 999 
CALL CONDI 
CALCONes 
CALICON DY 
Q= QO*ENP(-C2*(L:2-XNTOT)#*2)"EX P(-C2#(2.00 FW, 2 
& -TIMTOT#V)**2) 
EN1U.J)= TIME*(Q*DX + KI*T(I-1J)*DY (2#DXL) 
& +K3*TU+ LJ DY (22DNR)+ K4* TL. + 1)*DX DYB-(II*DX. 2 
& +KIPDY (2*DNL)+K3*DY (2*7DXR)+ K4*DX, DYB) 
& *TULJ)+ENUJ) 
See TOP WALL-NO FLUX 
ELSE IF (LEG.I) AND (ELT FLUXD OR iG he <2) 
& THEN 
999 SCALINGe In 
CALL .COnbe 
CALL Gon Da 
ENI(1J)= TIME*(DX*H*TI + K1*T(I-1.J)*DY/(2*DXL) 
& +K3*T(14+1,))*DY(2*DXR)+ K4* T(J + 1)*DX 'DYB-(DX*H 
& +KI*DY,(2*DXL)+ K3*DY (2#DXR)+K4*DX DYB) 
& *T(1J))+ENU,J) 
sess aes SE CORNER 
ELSE IF ((J.EQ.N+1).AND.(ILEQ.N+1)) THEN 
CALL COs DI 
CALINCON DE 
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ENM(1J)= TIME*((DX + DY)*H®TI 24+K1*T(I-1,J) 
& *DY (2*DXL)+ K2°T(IJ-1)*DX '(2*DYT)-((DX + DY) 
& *H 2+K1*DY (2*DXL)+ K2*DX ‘(2*DYT)) 
1 J)) + ENG) 
aeeanhes RIGHT WALL (NOT CORNERS) 
ELSE IF (I.EQ.N+1) THEN 
Gil COND! 
CALL COND? 
G1 COND 
ENI(I.J)= TIME*(DY*H*TI + K1I*T(I-1J)*DY DXL 
& +K2*T(1J-1)*DX (2*DYT)+ K4*T(L.J+ 1)*DX'(2*DYB) 
& -(DY*H+KI*DY DXL+ K2*DX‘(2*DYT)+ K4*DX (2*DYB)) 
& *T(I.J))+EN(IJ) 
aid BOTTOM WALL 
ELSE IF (J.EQ.N+1) THEN 
GAbbCOND|! 
Gal COND? 
GALL COND3 
EN((1J)= TIME*(DX*H*TI+ K1*T(I-1.J)*DY (2*DXL) 
& +K3*T(I+ 1J)*DY (22DXR)+ K2*T(1J-1)*DX DYT 
& -(DX*H+K1I"DY (2*DNL)+ K3*DY (2#DXR)+ K2*DX DYT) 
& *T(1Jn)+ENU.) 
feasts INTERIOR POINTS 
ELSE 
GALL CONDI 
GAM COND? 
GALIBGOND3 
CALL C@ND4 
ENI(1.J)= TIME*(K1*T(I-1.J)*DY, DXL+ K3*T(I+1,J) 
& *DY DXR+K2*T(L,J-1)*DX DYT+K4*T(I.J+1)*DX'‘DYB 
& -(K1*DY DXL+K3*DY DXR+K2*DX, DYT+ K4*DX, DYB) 
& *T(IJ))+ EN(LJ) 
END IF 
IF (£.GE.FLUX1).AND.(I.LT.FLUX2)) THEN 
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XTOT= XTOT+DX1 
ELSE 
XTOT= XTOT+DX2 
END IF 
205 CONTINUE 
200 CONTINUE 
FLAG=0 
DO 300 1=1,N+1 
DO 305J=1,N+1 
IF (1.GTELU XD ANDILT Flu) Tees 
DX = FW.(FLUX2-FLUX1) 
ELSE IF ((I-EQ.FLUX1).OR.(I.EQ.FLUX2)) THEN 
DX =(FW (FLUX2-FLUX1) +(L-FW)'(N + FLUX1-FLUX2)) 2 
ELSE 
DN= (EWN Fl ae) 
END IF 
IF (J.LT.NINT(.5*(N+ 1))) THEN 
DY =0.2"W/NINT(O.S%(N+ 1) 
ELSE IF (J.EQ.NINT(.5#(N + 1))) THEN 
DY =(.2*W'NINT(.52(N-+ 1))+.8#W (N-NINT(.54(N + 1)))) 2 
ELSE 
DY =0.8*W.(N-NINT(.5*(N + 1))) 
END IF 
C=418.485+0.2115*T(1J) 
ENMELT = TMELT*(418.458 + .2115*TMELT)*RHO#*DX*DY*S(I,J) 
IF (EN1(I,J).GT.ENMELT) THEN 
ENFUS = 247300*RHO*DX*DY 
IF (ENI(I,J).GT.SENMELT + ENFUS)) THEN 
T(1J)=(ENI(1,J)-ENFUS)(RHO#C#DX*DY*S(I,J)) 
K MELT(I,J)= FACTOR 
FLAG=FLAG+1 
ELSE 
T(1.J)=TMELT 
G K MELT(I,J)= 1.D0 


END IF 
ELSE 
T(1.J)= ENI(L,J) (RHO*C#DX*DY*S(L.J)) 
END IF 
IF (FLAG.EQ.0) THEN 
TIME=TIMORG 
ELSE 
TIME=TIMORG FACTOR 
END IF 
EN(1.J)=ENI(I.J) 
BOs (CON TINLE 
300 CONTINUE 
41S FORMAT(IN.26F5.0) 
C415 FORMAT(3F4.0,35F5.0,3F4.0) 
feo IF(CNI.EQ.1) GO TO'8P2 
IF (CNT.NE.INT(CNT 0010)*0010) GO TO 400 
C WRITE(1.*) 
$22 WRITE(*.*) CNT 
32] FORMAT(IN.FS.1.1N.F8.1,1N,F8.1,1X,F8.1.1X,F8.1) 
WRITE(1.421) TIMTOT 
DO 500 J=1.N+1 
WRITE(1.415) (T(1.J)-293,1 = 1,16) 
500 CONTINUE 
WRITE(1,*) 
C400 CONTINUE 
400 T91J)=TUJ) 
IF (TIMTOT.LT.15.6) GO TO 110 
eel) CONTINUE 
602 FORMAT(IX.1IAI1) 
CLOSE(1) 
C WRITE(*,*) QTOT 
END 
shubueee® SUBROUTINES FOR VALLES OF K 
SUBROUTINE CONDI 
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DIMENSION 1T(051,051),K MELT(051,051) 

REAL*S ALBLCEDILAG IT KIK2 RS hak MEET 

& DXLDXER Di Dis 

COMMON IIJLALBI-CIEDI, 1 Ki Ko Re Rae EL Ip xX 
& DXR,DYT,DYB 

INTEGER 

A=AI+Bl*T(1-LJ1)+ Cl*T(11-1,J1)**2+ DI*T(11-1,51)**3 


RI=A*BDXE- DAR (DNB DARA) ee rede 
IE S)B) 

SUBROUTINE COND? 

DIMENSION DOsh0sl) K VEL VOs10s1) 

REAL*S AL BIECTDIAB TAR VR2 RS Ra RAE De 

& DNR.DYT.DTB 

COMMOSR TIIUAPBLCRDI KIER? KR? ha ie ib. 
& DSRDr lpr 

IN FEGERSII SI 

A=Al+ BI Teer Jil} 92. Dias 


K2= A*B*(DYT+ DYB) (DYT*B + DYB*A)*K MELT(IL,J1) 
END 

SUBROU IIe Gor. 

DIMENSION 1(051,051),K MELT(051,051) 

REAL*8 AUPLChDIAB 1 KURO Raheh en, bar 
IDSC UDG Dn 

COMMON Th), AUBILGIDILT KIK2 co hoi ieee Dae 
& DNR,DYT,DYB 

INTEGER 11,1 

A=AI+Bl*T(11+ LJD+ClT(11 4 1,J1)**2+ DIA TI + 1J1)**3 
B=AI+Bl*T(11J1)+ Cl TL JL)**#2+ DIT LIL) 3 
K3=A*B*(DXL+ DXR) (DNL*A+ DXR“*B)*K MELT(IL,J1) 
EXD 

SUBROUTINE COND4 

DIMENSION 1(051,051),K MELT(051,051) 
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Reve os PIC DI-A.B,T.K1,Ke,K3.K4. KMELT,DXL, 

& DXR.DYT,DTB 

COMPO TONE B I Cl.D1. 1, Nee hok4 K MELT.DXL, 
oe DAR.DYT.DYB 

SiG Rott 


bP — eee OI Cl 1 Ji 2+ DIF (ll J1)**3 


K4=A*B*(DYT + DYB) (DYT*A+ DYB*B)*KMELT(I1,J1) 
END 
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APPENDIX D 
DATA 


The data is presented in the following format. 


Number on top line to the left is time (seconds) after start of data 


run. 


Data is presented as a temperature increase (K) over initial tem- 
perature (293K) 


Only the left half of the cross-section is presented; the right half 
is a mirror image for all runs presented. (The left column of data 
corresponds to the left edge of the cross-section. The right col- 
umn of the data is at the center line of the cross-section.) 


One-half of the vertical nodes are concentrated in the top 20 per- 


cent of the thickness and one-half of the horizontal nodes are 


concentrated in the first 0.01 m out from the center line. 


Each data run has four time sheets presented: 


1) 
2) 
3) 
4) 


initial heating phase 
initial onset of molten conditions 
maximum heat input 


end of run (torch completely past cross-section) 


Molten condition for aluminum: +641K 


Molten condition for steel: +1571K 
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Thickness / Current: .15m 170A 

Width / Welding Eff.: .o5m n=.8 
Voltage / Torch Speed: 28V v = 2.54 mm/sec 
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Thickness / Current: 15m 170A 

Width / Welding Eff.: .35m n=.8 
Voltage / Torch Speed: 28V v = 2.54 mm/sec 
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Width / Voltage: 35m 
Thickness / Welding Eff.: .15m 
Current / Torch Speed: 170A 
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Thickness: .15m 
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Voltage: 28V 
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Torch Speed: v = 2.54 mm/sec 
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Voltage: 28V 
Welding Eff.: n=.8 
Torch Speed: v = 2.54 mm/sec 
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Width / Voitage: 325m 28V 
Thickness / Welding Eff.: .15m Hn=.8 
Current / Torch Speed: | 153A v = 2.54 mm/sec 
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